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Abstract 

We return to the subject of stability of infinite time asymptotics of kinetic equations. We found 
a model which is simpler than those studied previously and which shows unstable behavior corre- 
sponding to our arguments from Q, where, however, a relatively complicated problem was treated. 
Our simplification to four levels interacting with surroundings enable us to proceed easily through 
all the way with just a pen and paper. We provide no numerical modelling whose justification 
causes naturally difficulties to the reader. We draw also further consequences of the found instabil- 
ity, not only with respect to higher order terms in kinetic equations but also concerning the very 
philosophy of physical modelling. The latter point can give more practically oriented physicist even 
better motivation than mere speculations about potential instabilities due to higher order terms 
in perturbation treatments without concrete resolution of correct asymptotics. 
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I. INTRODUCTION 



Subject of our paper are new results questioning stability of solution of kinetic equations. 
Kinetic equations are widely used in solid state or statistical physics in modelling of transfer 
processes including relaxation phenomena, influence of external fields etc. This concerns 
a great number of physical theories appropriate for different physical regimes of interest 
like Boltzmann equation 0, Fokker-Planck equation [[|, Pauli master equation |3J], or its 
generalized version introduced independently in different forms by Zwanzig, Mori |(J 
etc. General mathematical structure of these theories is the set of differential or integro- 
differential equations (of the first order), which determine time evolution of quantities of 
interest. Specific type of these differential equations is not unique, in its easiest form we 
meet with time independent Markov process (without memory), the other may contain time 
dependencies in coefficients (for example dependence on external fields), memory terms 
(time nonlocal equations), or inhomogenous terms or nonlinearity. Instabilities or chaotic 
behavior in case of complicated nonlinear equations are surely not very surprising. Below 
we will deal with the easiest form of such equations - the set of linear differential equations 
with constant coefficients. Specific topic of this work is scrutiny of stability of steady state, 
including that of the attractor nature of the steady state. Let us give some comments here, 
how this treatment is related to the other types of mathematical structure which are also 
used in physically similar treatments, and about physical consequences of this work. Time- 
local and memory-including theory of Markov processes are related in general (for general 
mathematical theory- see for example j7|) and also explicit equivalence for specific type of 
equations stemming from the Liouville equation could be found ||. The goal of this work 
is to argue that memory for Markov processes can be integrated into time local coefficients. 
In addition we are interested in steady state and infinite time asymptotics. Coefficients (in 
evolution matrix) often turn into time independent ones in infinite time region || [HJ. Some 
kind of Markov approximation here becomes exact. This provides the connection with our 
work. Influence of external fields is further property which can be often related with topic 
of our interest. 

The origin of kinetic equations is mostly in perturbational theories, which turn funda- 
mental microscopical physical laws into differential equations determining time evolution of 
macroscopic quantities available in experiments. For practical purposes the coefficients are 
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usually approximated by calculation of leading terms of the Taylor series in a perturbational 
parameter. Quite usual is using approximation of the second order in coupling to uncon- 
trolled degrees of freedom, what enables to incorporate connected relaxation phenomena. 

Then the time evolution of " quantities of physical interest" is solved. Formally, validity 
of all these approximations of the perturbational origin is limited on the time axis. On 
the other hand, these are the steady state and asymptotic limit which are of the greatest 
physical importance. Otherwise well applicable theory is thus asked to achieve reasonable 
results also for this time region. Expected results like the Boltzmann distribution were found 
in the simplest kinetic models, what gave physicists strong belief in general applicability of 
the particular kinetic theory. There is a statement that little changes of coefficients can 
not change in very dramatic way results of numerical studies, which is usually considered 
as "physically reasonable" . However, we want to argue that this statement is NOT true. 
Maybe some people (mainly mathematically oriented) can justify such a statement in general, 
but they mostly think, that it is the matter of purely mathematical toys that are seldom 
present in normal day physics. Our warning is: Nothing is so far from the truth. 

Recently we found in |IJ rather unexpected instability in asymptotical behavior in case 
of a not very complicated open system. The latter showed some not very standard features, 
but no direct indication of internal collapse in modelling of time development, of the type 
of, e.g., a clearly unphysical result. Many model approximations, where one can meet the 
below described instability, may seem to be quite usual or not very suspicious. 

However, whatever is the truth about correct time asymptotics of the model hamiltonian 
from UU, we would like to warn the reader against unexpected features encountered in 
this problem. When one incorporates (according to his/her opinion) important physical 
processes into kinetic equations, attention must be also focused on stability of the solution 
of the problem. In particular the question is which part of the results is correct and which is 
only a belief (and contingently only an unjustified belief). In case when the resulted steady 
state can not be confronted with, e.g., thermodynamical laws, the result of thus invalid 
simulation may seem to be quite good. 

Sometimes physical belief in the second order approximation results is hidden into so- 
phisticated mathematical methods which pretend to give a "general" proof of correctness 
of the second order results without realistic treatment of the stability, case by case. Such 
methods might, regardless of their mathematical validity and usefulness, in concrete appli- 
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cations leave a space for speculations about existence of instabilities, in some important 
features of the solution, like in our model from [|TJ. In []J we showed how it is possible to 



fulfil so called Davies theorems [IT] and, simultaneously, also obtain sharp instability in the 
long time regime behavior against higher order contributions to Tokuyama-Mori coefficients. 
The proper quantity of our interest is time asymptotics of the kinetic equation (regardless 
of whether the steady state is given only by internal dynamics, or it is influenced by some 
regular external field, presence of interchange particle etc.). 

Organization of our paper is the following: In part 2 we introduce a quite usual transfer 
process model - four site open quantum system interacting with surroundings. We perform 
standard treatment using second order perturbational theory of relaxation with full spectral 
analysis of evolution matrix and the unique asymptotic steady state will be found. One can 
verify that the model, at least in standard thinking, contains description of all the important 
physical processes present here. Further we introduce (in section 3) a correction to evolution 
matrix. We thus add a transfer channel, whose description is constructed in the same way 
like that of the previous ones. However, we assume the effectiveness of the new channel to 
be very small with respect to that of the original channels. This new process is also not 
apparently seen to change overall properties of system. Nevertheless, dramatic changes in 
the steady state will appear. The small perturbation will be treated in two distinct ways. 
First we consider it like a higher order correction in the perturbational series to our second 
order approximation. This approach is in our view mathematically more interesting, but is 
rather speculative because we do not explicitly calculate all the fourth order contribution 
that can come from higher order calculation in some exact microscopical theory. However, 
the instability of the second order approximation is verified regardless this objection. The 
second possibility is to complement the model by a further small parameter which will be 
limited to zero after the asymptotic calculation. In this case we fully stay at a position of the 
standard second order kinetic theory, but we consider model character of our hamiltonian 
as in the everyday physics, and we treat stability of the model conclusions against "omitted 
physically irrelevant" process. In part 4 we throw some light on indication of instability of 
the given type. At first we show its indication in spectrum of the evolution matrix, then 
a purely analytical treatment of stability will be given. In part 5 we compare our result 
with van Hove limit of the same model Hamiltonian. We will find a connection with some 
pathological case of this well understood limit, where kinetic theory also may have problems 
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with establishing of asymptotical state. 

All these ideas are presented on this simple model. Further consequences about relation to 
Davies theorems are not repeated , interested reader can find them in [|IJ . All the calculations 
are hand made product and the reader is pretty invited to follow them. The reader is 
also invited to think how the referred problems, which one may consider as clear or quite 
trivial, may become forgotten when a complicated model is treated which is not analytically 
calculable, but where only computer simulations are at hand. 

II. MODEL 

Let us consider the following model hamiltonian of a four site open system. We measure 
energy in units K. 

H = e(c\ci+clc 3 )+J (4 c 3+4 c 2)+XX ( ^i 1 ~^ (- B fc c i c 2+-Bfe4 c i)+£fe 3 ~ 4) (B k (^c±+B\(^cz)+VL k B\B k .} 

k 

(1) 

This hamiltonian describes our four site system where creation cf and annihilation operators 
Ci are each related to the i-th site, with three transfer channels, two of them being bath- 
induced the remaining one being coherent. As far as we consider one particle only, there is 
no necessity to introduce (anti)commutational relations between these operators. Dynamics 
of the system is dominated by bath induced transfer channels between 1-2 and 3-4 sites and 
a coherent transfer channel 2-3. B k , B k are creation and annihilation operators of k-th bath 
phonon mode ( fulfilling boson commutational relation), G k are coupling constants of the 
system-bath interaction. Parameter J describes power of the coherent channel. 

We do not allow the interference between bath induced channels 1-2 and 3-4, what means 
to fulfil conditions like: 

]T 6(e - n k )Gt 2> Gt 4) Tr Bath (p Bath BtB k ) = 0. 

k 

This can be generally fulfilled if one considers that the particular transfer channels are 
induced by different phonon modes 

^(1-2)^(3-4) _ n 

One can treat this hamiltonian using various schemes. Firstly, it is possible to think 
about different regimes, according to different magnitude of the coefficients in Hamiltonian. 
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We are here interested in the regime where it is appropriate to treat J- and bath-induced 
transfer channel as a perturbation. We emphasize that this choice does not correspond to 
so-called van Hove limit [|12| . One can diversify the physical interpretation of hamiltonian 
(|I|) and a chosen perturbation scheme. Coherent channel one can consider as a slow internal 
motion treated according to [IB], but it may also represent constant or periodical external 
field influence. 

Various constructions of kinetic equations can be also applied. We restrict ourselves to 
those which respect chosen mathematical structure and the physical regime. Though also 
here physicist use various formalisms, one may obtain our results using Nakajima-Zwanzig 
identity, Tokuyama-Mori equation (both in their second order approximation), and also 



Haken-Strobl-Reineker parameterization |14|, [15]; all these ways lead to formally the same 
master equation: 



dt 



J2 w {ij},{ki}P{ki} 



(2) 



{hi} 



where vector p is organized in the following way 

P T = (pn, P22, P33, P44, Rep2z, Imp 23 , Rep 12 , Impn, Repis, 



Imp 13 , Rep M , Imp u , Rep M , Imp 24 , Rep u , Imp u ) 
Matrix W we call evolution matrix. is the second order approximation of W. It reads 

/ A n n n \ 



A 














B 














c 














D 



(3) 



where 



.4 





r t 
















-r T 











-2J 








-ri 


r T 





2J 










-r T 




















2 


— e 





J 


-J 





e 


9 



(4) 



6 



B 



C 



2 

— e 

J 

2 

— e 


-J 



6 

r T +ri 

2 

-J 





e 

_£i±£i 

2 
J 



r T +r x 





J 

-r ; o 
o -r 



j 



(5) 



I J 



D 



2 

— e 





-J 

_r T 
o 

e 

r T +ri 



J 



-r 



(6) 



(7) 



Here 



r T = 2TxY J [G { t 2) f^-^k)Tr BathPBath {BlB k ) = 2nJ2[G i t i) } 2 ^-^k)Tr Bath (p Bath BlB 



(3-4) 1 2; 



= 27r^[G , ^ 2) ] 2 5( e -fi fc )Tr i?a ^(p i?a ^ J B fc 4) = 27r^[G , r 4) ] 2 ( 5(e-fi fc )Tr i3ai ,(p i3ai , J B fc 4). 

(8) 

The equality of coefficients for 1-2 and 3-4 transfer is our additional assumption, that can 
not be deduced from (p. Notice that J, T^Tj we consider as perturbations of the same 
magnitude, proportional to the parameter A 2 . 



J, T| oc 



(9) 



(Reason the proportionality is only a consistency with standard perturbational order of 
bath- induced transfer channel.) 

The steady state is given by the condition: 



W {ij},{kl}P{kl} = 

{fci} 



(10) 



We are now to calculate the complete spectrum of the evolution matrix. Firstly: this 
enables us to show that steady state is also the unique asymptotic state of this equation. 
Furthermore we will argue that the solution has no apparent deviant feature. Last but not 
least, in section 4 we will show that in a careful treatment one can indicate, in this spectrum, 
the instability calculated below. 
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The evolution matrix was arranged so that it has a quasidiagonal structure. We have to 
calculate a characteristic equation. After bit of algebra (we must calculate determinant of 
submatrices of maximal order 6) and rearranging resulting terms, we obtain 1 : 

o = £ • + r T + r ; ) • m + r T + r ; ) [(£ + + e 2 ] + u\t + T -^^) 2 } 

■ft + ie+ + Ti) + J 2 ) ^t^XZ + ?i) + J 2 } 

+ ^ + — + Ft) + j2} ' { ^ ~ ie + + Ft) + j2} 

•Ke + ^^^ + e 2 ] (ii) 

Twelve roots can be calculated directly from the quadratic terms. What remains is an equa- 
tion of the fourth order. The roots can be in principe also extracted using Cardano formula, 
but it does not provide an easy survey. Instead we inspect behavior in the A — > limit of 
the perturbational parameter. This analysis and calculation of twelve exact eigenvectors is 
provided in Appendix. 

In conclusion: there is only one steady state and, at least for not very high parameter 
A, all the other eigenvectors of matrix |3| have negative real parts, i.e. connected terms 
in time evolution simulation disappear in the infinite time and the steady state is also 
the asymptotical one. Because of finite order of the matrix there is a region surrounding 
the zero, where this problem has infinite time asymptotics given by ([TOD, so one can limit 
himself /herself to this region without complications. For very high values of A parameter, 
the model need not have the correct behavior in accordance with general inapplicability of 
the perturbational treatment for this case. The worth notice is that the zero eigenvalue 
came purely from the first submatrix A. The others have nonzero determinants and thus 
the only solution of the steady state condition must be zero for associated elements in the 
density matrix. 

We can give at this place the solution to steady state condition fllOD. We work with 
normalization condition 

£Pn = l- (12) 



The terms are ordered according to ordering of submatrices; one submatrix is one row. 
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Then the result is: 

pn ~ (r T + r^ ; P22 " P33 " (r T + r^ ; P44 " (rVTT^ ; (13) 

We specifically note the equality in population at sites 2 and 3. The reader may speculate 
whether this model and result are in whatever sense bad. In any case, there is no internal 
collapse in these calculations. One may be suspicious about the fact that the model does 
not lead to thermodynamical equilibrium, but as we notice above we do not argue that this 
system is in thermodynamical equilibrium as we are not in the van Hove limit, and also 
the purely internal character of the coherent transfer was not specified. There are some 
physicists who believe that the validity of thermodynamical laws must be in some direction 



connected with the van Hove limit [16|. The 2-3 channel is elastic what implies 2-3 symmetry 
and consequent equality P22 = P33- The transfer term proportional to J can come also from 
defined (e.g. harmonic) external field; then thermodynamical prediction can fail or this 
prediction need not be clear without further calculations. 



III. PERTURBATION 

In this section we introduce a small perturbation of model (fj) in form of an incoherent 
transfer channel between sites 2 and 3. The new terms in evolution matrix can be quite 
small with respect to the other ones coming from the previous consideration. Construction 
of the terms is fully analogical to the previous one. Formally one can include a change into 
hamiltonian : 

5H = Y,Gt 3 \B k cic 2 + Bi4c 3 ). 

k 

We refer a new corrected evolution matrix, the effectiveness of the new 2-3 channel is mea- 
sured by rate coefficients g^Qi- 
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Now we calculate the stationary state of problem (|T6|) . We again work with the assumption 
that A is so small that submatrices B, C and D are regular and so the stationary condition 
applied here has the trivial solution only Then one gets the result after an easy algebra: 



P22 = CT l T^{g l 
P33 = CT^T^{g^ 
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where C is a normalization constant to be deduced from (O): 
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We are interested especially in the ratio of the population on the sites 2 and 3. The 
reason for this specific interest becomes apparent later. 
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(17) 



The term "small perturbation" has to be formalized in order to talk about the instability. 
We have worked out this point in two different ways. First of them is submitted mainly for 
a mathematically oriented reader. We consider g^,g[ as proportional to A 4 . 

j,r i5 r T ocA 2 ^ T ,^ocA 4 



(18) 
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One may have some physical objections against this interpretation, stemming from the fact, 
that we did not provide complete 4-th order inspection of the kinetic theory. But we have 
quite narrow ambition here. We point out the instability of the result ( pT3[ ) against the 4-th 
order correction, which we consider to be arbitrary - as a potentiality. We argue that an 
arbitrary perturbation could be used to achieve this conclusion. The physical motivation 
in this interpretation is in the background only, in order to get the reader interested, the 
statement is of mathematical character. One can also omit here the additional term in 
Hamiltonian 5H , and think of the perturbation as of a higher order terms coming potentially 
from the Hamiltonian ([!]) that are omitted in standard second order calculation. We will 
compare the results (|13|) and (|l6l) in A — > limit where the perturbational treatment is best 
verified. (Performing this limit has of course no consequence in connection with the main 
statement - instability.) 

The second interpretation stays fully on the position of the second order kinetic equation. 
We introduce some further parameter rj that measures relative power of different transfer 
channel 

^r h r T ocA 2 g^ocr/A 2 . (19) 

After evaluation of the A — > limit, which gives precise mathematical sense to our calcula- 
tion, we consider rj to be small, formally limiting to 0. We shall show that regardless the 
arbitrarily small (but nonzero) magnitude of rj, the result ([13]) is not preserved. In other 
words, 

limp(?7) p{rj = 0) 

rj— >0 

p designates here the steady state in A —>■ limit. This is the central statement that we 
are going to prove. We inspect the stability of standard kinetic equations with respect to 
physical processes which were not incorporated into a model in question because of their not 
high strength (at least from a formal, cursory point of view) and consequent underestimating 
their influence. This point is possibly not so interesting mathematically because the second 
order theory is held here, but it seriously questions the straightforward applicability of 
the standard kinetics from the physical point of view. Both these interpretations are from 
mathematical and also from physical context quite distinct. We argue that the instable 
behavior is the internal problem of the approximation (|^) and does not come from the very 
specialized choice of perturbation or scheme of its treatment. In the next subsection we 



11 



make this point more clear. 
Performing the announced limits: 
First interpretation (|18|): 

\4_ I 2£ , A^rt+r^+A^gt+g^ N , 9X 4 72 

l im PH = i im ^ u A2(r T +r i )+A4(g T +g i) + 2 ) + AA J = g± 

-0p33 ^0 A 4, t( v(rT+r + A'CTr^HA^) ) + 2A4J 2 * 



Second interpretation (0): 



ri \2_ / 2e 2 1 A 2 [r T +rj+??(gT+9|)h , r>\4 72 

lim lim ^ = lim lim ^ ^^^^ A2[r +r ! f . J = lim ^ = ^ 

^OA-,o P33 ^ ^ TyA^ T ( Aa[rT+r ^ faT+gi)] + £W±toill ) +2A^ *-°ffr St 

Also other results are identical in both our limits, we refer it in a short way. 
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.Re^ = 0; Imp 23 = (20) 



We clearly see that if ^ g±, then sharp changes appear in solution (|IBD. We did 
not assume equality between these coefficients. Moreover the physical motivation of our 
correction points out to ^ g\. Rather the standard thermodynamics relation 

m n 

— = exp pe 

9l 

may be assumed. Derivation of this statement consists in some additional assumptions about 
initial state of the bath, which are, however, standard. Of course, if one is interested in the 
Taylor series structure in higher order expansion of (|T7|) , our interpretations are mutually 
different, but the general picture is not changed. 

We would like to give further warning here. One need not be very surprised because of the 
following argumentation. Transfer channel connected with parameter J is in fact physically 
also of the fourth order, because its direct application to density matrix (commutator [J, p\ ) 
change either the ket or bra side of density matrix only and comparable process connecting 
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the diagonal terms in density matrix is thus of the fourth order. Then this new included 
channel is comparably strong (in the first interpretation) or even stronger (in the second 
argumentation) than the first one. We give threefold counterargument: 

1, - Nevertheless, formally the coherent process (J-proportional) is included in the second 

order. Such a treatment is absolutely standard. Care in this direction becomes difficult 
or technically unable for complicated system. Moreover this objection also seriously 
questions concept of the generalized master equation in general, because treatment of 
the whole density matrix (of the system) included information about set of generally 
incompatible observables. With respect to a particular measurement (here for example 
site probability measurements) the other terms unrelated to this measurement (here 
off-diagonal matrix element) have always the role of some kind of memory. So their 
treatment apparently differs from that of bath induced channel one. 

2, - Physical intuition in more complicated case is uncertain and may fail. 

3, - If one is internally sure about his/her intuition, please try inspect the formula (|17D 

again now with the following sheme 

T h T l ,J<x\ 2 ; #1,^ oc A 5 

or alternatively 

r T ,r b JocA 2 g h g i oc r/A 4 . 

One can see here that in the same limits as above, the instability is still present, though 
the perturbation should be now smaller than the fourth order coherent channel. In 
fact the coherent channel transfer is seemingly of the higher order then 4! Have you 
seen it before this limit calculation? 

IV. INDICATION OF INSTABILITY 

We have argued in the previous section, that instability described above is the matter of 
internal problem of approximation (|3]). Despite of its reasonable behavior in time evolution 
simulation which we proved in Appendix, a problematic step has been indicated above and 
this problem must be visible purely from the second order approximation calculation. In 
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fact, if one were not able to give some indication of the instability of the approximation 
from itself, the situation would have been critical. At first: Calculation in higher order 
is not standard, and what is worse, it is an extremely difficult task. The second problem 
is that in that case one can never ( in no finite order of calculation) be sure whether the 
provided approximation is already stable. We give two points which are connected with our 
instability and indicate it. 

A. Long-time excitation in the spectrum 

This method enables us to show quite simply how to obtain the indication purely from 
spectrum of the evolution matrix in the second order approximation. What is important 
from practical point of view, one may use this method without principal difficulties also in 
complicated models by numerical analysis and, further, may in many ways also implement 
orientational indication into her/his time evolution computer simulation. Only for simplicity 
and without any change in physical context we assume that the number of linearly indepen- 
dent eigenvectors equals to the order of the matrix. In problem this statement is proven, 
because none of the submatrix has two identical eigenvalue. The general form of spectral 
decomposition of finite (non-normal) matrix is then: 

q 

where £ g is q-th eigenvalue and R q (L q ) is the associated right (left i.e usual) row (column) 
eigenvector. This decomposition holds good this normalization: 

i 

Kinetic equations conserve total probability, what results into the fact that eigenvalue 
is always present. Have a look at our result in the Appendix again. The very suspicious 
eigenvalue is £4. In perturbation scheme (§), there is a proportionality £ 4 oc A 6 . No wonder 
that this eigenvalue need not be very stable against our perturbation regardless of the type of 
scheme. More generally: Let us have, in our spectrum, eigenvalue with real part approaching 
zero (with A — > 0) and proportional to higher than second power of A (n > 2): 

W tj = • (L ) i (R ) J + ai X n ■ {L x ) i {R 1 ) j + ... 
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We explicitly emphasize that also (L g ),(R q ) are A-dependent what only enables that 
has terms of just the second order. 

Then one can easily construct the mathematical "perturbation" which causes the insta- 
bility, for example perturbation 

SWij = -X n ■ (LoUR^j - ai X n ■ {Li)i{Rx)i + 0(X n ) 

that change the stationary state from (Lq) to (Li) in the limit. 

In the first perturbational technique fll8l) possibility of such construction is straightfor- 
ward. In the second case ( |l9l) we must construct it from the more complicated terms, but 
one can see in our model ( |l6l) that it is possible. Of course, not each one of the math- 
ematical "perturbations" is physically interpretable. One usually has some conditions for 
the evolution matrix stemming from the conservation laws (at least particle conservation in 
case of solid state physics), etc., but the set of possible perturbations is so great, that it 
surely contains also reasonable perturbations. One of them we introduced in subsection 3. 
In our problem, there is the " near-the-zero" eigenvalue proportional to the sixth order in A, 
so we could choose the perturbation smaller than we have done. Concluding this subsection 
we notice that the reciprocal real part of eigenvalue (without sign) can be also called the 
lifetime of "excitation" (though it need not be a very appropriate name in some cases) . 
This provides practical indication of this instability - highly increasing lifetime of relaxation 
phenomena when the parameter of perturbation is reduced according to formal scheme of 
construction of given kinetic equation. Because of clarity of this point we do not give any 
formal statement. 

B. Analytical treatment of stationary condition 

In this subsection we give an easy example of analytical stationary condition flT0| ) treat- 
ment. This treatment is stable against perturbation. We know that it need not be usually 
the very appropriate method from practical point of view. When a direct explicit resolution 
is not available (for more complicated or extended problems) one must take extreme care in 
computational implementation about numerical errors. The main reason for introducing this 
calculation is further understanding of the origin of the instability for the reader who still 
has not internally accepted the presented facts. For our treatment we need (for simplicity) 
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to assume that there is a unique asymptotical state of the system. This is because we will 
here take care of the stability of the zero eigenvector of the evolution matrix only. This 
treatment does not provide the proof that there is no potential eigenvalue with positive real 
part (collapse of model) or there is some eigenvalue along the imaginary axes so near to zero 
(real part) that it can approach through some perturbation the imaginary axis. We look for 
resolution of approximation (|3J) as Taylor series coefficients. 

p = J2 \ n p {n) 

n 

The important difference as compared to Taylor series of solution (|13j) is that we explicitly 
assume existence of perturbation of the order A 3 , respectively 77A 2 , which is otherwise arbi- 
trary. We take only the results which are independent of potential perturbation. However, 
this is a standard correct perturbational method. Such a treatment gives us only finite 
number of conditions for Taylor coefficients. The calculation is straightforward. Condition 
in the zeroth order enables the calculation of 

p!?=0; pff=0; $ = 0; p£? = 0, 

while in the second order 

Jo) _ n - n (°) - n- - n- - n- t™ - n- r p J 2 ) - t(J® n^h 

Pl3 - U > P24 - U > Pl2 - U ) Pu - U ; lm P23 - U > ^ e P23 - J {P33 ~ P22 ) , 

= P ( 2 = ^/®. (21) 

Here we clearly see the internal problem of the second order approximation (|3|). Neither 
the zeroth order of the density matrix is resolved by stationary condition (|T0|). We have 
still a two-dimensional subspace (arbitrary p^'iPm) where the steady state can be found 
(in the Liouville space of the density matrix). The result ( [TBI ) is the corollary of implicit 
assumption of zero effect of higher order calculation, not justifiable from the mathematical 
point of view. One can comprehend that including a potentially higher order perturbation 
like (p~6|) will define the zeroth order density matrix in space of our result ( pT[ ) with a high 
degree of arbitrariness. 

The solution, which is correct to at least the zeroth order and would have to be stable 
against perturbations, one must calculate more precisely. 

We notice that including the perturbation not as a potentiality, but like a real effect 
changes the stability of the model. One can prove in both interpretations that we obtain 
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further condition (once, in the order A 4 or respectively r]X 2 ) 

(o) (o) 

#tP22 = 91P33 ■ 

In the first interpretation we obtain in the fourth order also further conditions; however, 
because of its speculative character we will not publish it here. This means that the situation, 
despite of its unpleasant character, is not hopeless. One can indicate instability and also 
the ways to improve models are principally possible. Let us notice that the result obtained 
in this subsection is in accordance with all previous calculations in A — > limit. 



V. THE VAN HOVE LIMIT 



All our previous results were obtained in a way that is not just standard in relaxation 
theory. We introduced perturbational scheme (|9|) for calculation of the second order kinetic 
equation for model Hamiltonian ([TJ). We gave some physical arguments for this choice. 
Nevertheless, the standard variant of great popularity is of course the van Hove limit ||12||: 



J oc 1, T T , T l oc A 2 . (22) 

We argue here that the problem with infinite time asymptotics of the model ([!]) (in the 
second order kinetic equation) is reflected also in this well understood limit. To see this we 
introduce " energetic" representation in eigenvectors of H$ 

c„ = ac 2 - pc 3 , c m = ac 3 + /3c 2 (23) 

where: 

vgl + y/r+^) V2J J 

a = — i = oc 1, p — — -, = oc — 

2\ hS- + 1 + v/ifTI eM + 1 + WTi e 



The model (|l|) is now: 



e / 4J 2 e / 4J 2 

H = ecld + -(1 + Jl + _) c { n c ra + -(1 - \jl + —)c\ I c ll + Y / {n k BiB k + 



' *- k 

G^- 2 \B k c\{ac ll +f3c lll )+Bl{ a c\ l +(3c\ l ^ 

(24) 



We offer some comments regarding Hamiltonian (24) and the van Hove limit (0). There 



are two bath induced channels between levels 1,11 and 111,4 respectively, in analogy with the 
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previous treatment in site representation. What is the difference is that there is no coherent 
transfer term in energetic representation; on the other hand two weak bath induced channels 
between levels 1,111 and 11,4 appeared. The strength of these channels is proportional to 
(JG) 2 , so this term is in the second order kinetic theory of relevance in the van Hove limit 
only. (The region of physical applicability of does not contain the regime specified 
before in connection with (|9|). We only clarify behavior of treated model hamiltonian from 
another view.) These channels cause communication between specified levels for short time 
regime, nevertheless both the channels lie off the energy shell, so for the long time regime 
this transfer is forbidden. Then the second order kinetic theory with integrated memory 
like in infinite time forbids all the communication between pair of levels l+II and that of 
levels III+4. Asymptotical stationary condition then has two linearly independent solutions. 
Nevertheless 

c\ci + cjjCn 

does not commute with the full Hamiltonian (|24]), it is no integral of motion. Consequently 
one can not use the long time (Born-Markov) approximation upon looking for time asymp- 
totics - the result may depend also on short time transient effects. The result obtained in 
this way is also seemingly unstable against higher order calculation. In Appendix B we 
give the complete second order kinetic equation and its solution in the van Hove scheme. 
The solution ( |B"6| ) of stationary condition (|10j) shows just the same asymptotic state of the 
density matrix (and potential instability) like fl2T|). 

We stay here at quite real physical problem: Consideration whether these levels are 
isolated and the transfer is strictly forbidden, what suggests ordinary meaning of the energy 
conservation law, or whether a limited value of electron density can be transferred. The 
significance of the van Hove limit is here also questioned. 

Last but not least, we may also think about more symmetrical case of system - bath 
coupling of form: 

J2{G { t 2 \B k + Bt)(c\c 2 + 4d) + Gf\B k + 4)(4c 4 + 4c 3 ). 

k 

The result for ([|) remains unchanged. Then in energy representation also 1,111 and 11,4 
channels appear that lie on energy shell for interaction with low energy phonons, if these are 
present. Then the asymptotics in the van Hove limit can be obtained as unique. However, 
the existence of such phonons (and channels) is a serious change in physical meaning of the 
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entering model. Inapplicability of the van Hove approach for physical regime J < G 2 is 
apparent. 

VI. PHYSICAL CONSEQUENCES AND CONCLUSIONS 

Let us firstly discuss some physical consequences of possible instability in the first scheme 
(|18|). We have shown that the result obtained in the second order may be unstable against 
some correction arising in a higher order calculation. Let us stay on the position that 
evolution generator (Hamiltonian) (|l]) is exact, and the only question is the correctness of 
its study. Of course, in that case, our construction of higher order correction is rather 
speculative, but it proved the instability of the result. The full analysis of the spectrum 
of the transfer matrix showed that the most 'slow excitation of the steady state' calculated 
in the second order has its lifetime of the sixth order in A. This gives that also in the 
case where first correction to terms W{22}{33}, W{33}{22} is of the order six, it may cause 
such an instability. This questions some standard results obtained in everyday simulations. 
Unfortunately we have also some intuitive physical arguments that one must really calculate 
at least six order processes for achievement reliable result for transfer between 2-3 sites. This 
is because the transfer rate really IS a process (at least) of the 6th order in specified scheme 
(§). Incoherent transfer term between sites 2 and 3 lies off energy-shell, consequently some 
collaboration with bath modes for stabilization is necessary. All the Feynman graphs one 
can draw must have for such a transfer great number of lines in order to get on energetic 
sphere, implying high order of this transfer. The reported delicate situation in the van Hove 
limit also supported caution against straightforward use of kinetic theory. Unfortunately, 
any higher order contribution calculation is a difficult task, getting dramatically worse from 
order to order. This fact will also in future cause the great popularity of the " naive " 
treatment; in this direction, our expectation concerning influence of our work is rather 
pessimistic. Greater care about applicability of usual model methods is, on the other hand 
and in the light of our results, more than appropriate. Further investigation should be turned 
to higher order inspection of master equation connected with Hamiltonian ([[]) in specified 
physical regime. Especially the question connected with the infinite time asymptotics is a 
great challenge, of crucial physical implication, and not satisfactorily resolved yet. 
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APPENDIX A: DETAILED ANALYSIS OF SPECTRUM 

Determination of the twelve eigenvalues does not meet with problems, because these are 
the roots of quadratic polynoms. In addition we would like mainly to know signs of real 
parts of the eigenvalues, at least in the limit A — > 0. Thus we reduce complicated results 
into the Taylor series at least to the order which gives the sign. : 



6 = 0, & = -r r 



.e 3r| + r T . e 2 . r,-r T r,-r T 2 T9 
f 7 = -i + i\ — + it— - — — + j ~ -r, 

s 2 4 V 4 4 4 



• e 3r ; + r T . /e 2 ,.r ; -r T (r\-r T ) 2 T2 . r ; + r T 



• e 3r T + r, ./e 2 . r 4 -r T (r^-r^) 2 T2 . r 1 + r l 



e + . e 2 . r^-r^ (r^-r^) 2 T2 

Further eigenvalues are roots of the 4-th order polynomial coming from the submatrix A. 
Though there is a formula which enables explicitly to extract the roots - so called Cardano 
formula, we do not use it because of its complicated form, and we only determine leading 
terms of the limit case A — > using the Taylor series. (This point provides no additional 
assumption about analytical structure of this dependence, all the results can be proved using 
mean value theorem.) 

f ~ r r- p ~zZ!iEi±Ii2 

?3 ~ -1 T ~ 1 h 44 ~ ~ 2 
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rV+Ti. ,* 

o > « - ?5 



Notice: The complex square root used in formulae above is defined into the upper half- 
plane of the complex plane (e.g. Imy > ). 

APPENDIX B: TIME ASYMPTOTICAL SOLUTION OF THE SECOND ORDER 
KINETIC EQUATION OF MODEL IN VAN HOVE LIMIT 

We start from (p^ ) and in the van Hove perturbational scheme (|22]). Organization of 
column vector of the density matrix is following: 



P 



Pu, Piiji, Piiijii, P44, Rep ni n, Impiijn, Repiji, Imp XIh Rep 1IH , 



Impijn, Repine Irnp niA , Rep HA , Imp IIA , Rep 14 , Imp 14 
Kinetic equations (0) obtained here from e.g ( 



where 



,_J are: 

f A ^ 

5 

C 

D J 



A 



pt; 

1 1 


pu 
1 T 








6T^ 





pu 

1 1 


pi) 
1 T 



















pi) 


pt; 

1 1 


-erf 










0F2 


pi' 

1 1 

*rj 


pu 

-1 T 


T + 1 




— e — 


2 


2 


2 


2 


2 



2A 



2 

-e - A 
0r| 





e + A 



T I 



2A 



2 w 

o S 

2 







p« 



-A 



a -ry, 



(Bl) 



(B2) 



(B3) 
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where : 



C 



2 

-e - A 



e + A 

pn I pi; 



9 Ii \ 

2 u 



u 2 



2 





-1 T 



-A 



2 

\ - 6 



pi) I pi) 

1 T + J- 







/3 



4J 2 
2J 



1) 



4J 2 ' 



= 2na 2 J2[G { t 2) } 2 ^ + A - n k )Tr Bath p Bath (BlB k ) 

k 

= 2WJ2[G k 3 ~ 4) ] 2 5(e + A - n k )Tr Bath (p Bath BlB h ) 

k 

r I = 2 ™ 2 E[ G ^ 2) ]^( e + A - n k )Tr Bath (p Bath B k Bt) 

k 

= 2ira 2 E[Cr 4) ]^(e + A - fi fc )Tr Bat ,(p Bat ,5 fc 4). 



(B4) 



(B5) 



One can verify that stationary condition (1C ) is satisfied by density matrix: 

~pv ~pv pu pi; 

P = ^(f^fp^ 1 + ^TT^ c n Cn ) + ^ " ^)(^fp^ c " lCin + p^Tr^ (B6) 

with arbitrarily chosen constant C G (0, 1). 
This proves the statement of main text. 
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